AD-A231  049 


AHn  azlTSSi-.ih/h/) 


% 


CONTACT-IMPACT  ALGORITHMS 
FOR  PENETRATION  STUDIES 


Ted  Belytschko 
Wing  Kam  Liu 
Mark  O.  Neal 


Nov.  1990 


U.S.ARMY  RESEARCH  OEHCE 
DAAL03>87-K-(X)35 
NORTHWESTERN  UNIVERSITY 


ELECT c 

JAN  18 1991 


APPROVED  FOR  PUBLIC  RELEASE; 
DISTRIBUTION  UNLIMITED. 


UNCLASSI7IED 

iS«bA»r?7'13:5Sl?)'OT<5'^"(5r’TK»«  '3a(5S 


>1ASTER  COPY 


FOR  REPRODUCTION  PURPOSES 


KEPOrr  DOCUMENTATION  PACE 


1«.  SEPORT  SECURITY  CLASSIFICATION 

- - - 


RESTRICTIVE  MARKINGS 


2*.  SECURITY  CLASSIFICATION  AUTHORITY 

2b.  OEOASSiFtCATION/OOWNGRAOINC  SCHEDULE 


3.  OlSTRiaUTIQN/AVAiLAglUTY  OF  REPORT 

Approved  for  public  release; 
distribution  unlimited. 


4.  PERPORMINC  ORGANIZAHON  REPORT  NUMSERCS) 


S.  IMONITORJNG  ORGANlZAHON  REPORT  NUMBERIS} 


6a.  NAME  OF  RERFORMINQ  ORGANIZATION 

Northwestern  University 

6b.  OFFICE  SYMBOL 
(tf  tpptkabtai 

7a.  NAME  OF  MONITORING  ORGANIZATION 

U.  S.  Army  Research  Office 

1  6e.  ADDRESS  (Ctf,  Stan,  and 

Civil  Engineering 

Evanston,  IL.  60208 

7b.  ADDRESS  (Oty,  Stan,  and  DPCodt) 

?.  0.  Box  12211 

Research  Triangle  Parle,  NC  27709-2211 

8*.  NAME  OF  FUNDING/ SPONSORING 
ORGANIZATION 

U.  S.  Army  Research  Office 

9b.  OFFKE  SYMBOL 
Ot^bpMcabtai 

3c  ADDRESS  (Ctr,  Stan,  and  DP  Csdbi 

10.  SOURCE  OF  FUNDING  NUMBERS 

P.  0.  Box  12211 


PROGRAM 

PROJECT 

TASK 

WORK  UNIT 

NC  27709-2211 

EUMENT  NO. 

NO. 

NO. 

ACCESSION  NO 

Contact  -  Impact  Algorithms  for  Penetration  Studies 


12.  PERSONAL  AUTHORfS)  ^  ^ 

Ted  Belytschko, 

Wing  Kam  Liu,  Mark  O'Neal 

I3a.  TYPE  OF  REPORT 

FINAL 

13b.  TIMS  C 
FROM  3 

22 

8/90 

15.  PAGE  COUNT 

76 

16.  SUPPUMENTARY  NOTATION 


17.  COSATl  CODES  I 

FIELD 

GROUP 

SUBGROUP 

The  vlevy  op In lone  and/or  findings  contained  in  this  report  are  those 
consc^ed  as,  an  official  Department  of  the  Armv  oositior 

orHo-r  r  ^ 


It.  SUBJECT  TERMS  (CoiTCinue  on  rw¥orm  if  nvctaury  jntf  idonatf  bf  biedt  nwnoor) 

Impact,  Penetration,  Finite  Elements 


ASSTRACr  (Condnue  on  rw¥orm  it  neonury  tnd  idonbff  by  biode  mmbor) 


Contactr impact  algorithms,  which  are  sometimes  called  slideline  algorithms,  are  a 
computationally  time-consuming  part  of  many  explicit  simulations  of  nonlinear 
problems.  One  of  the  reasons  for  their  computational  expense  is  that  slideline 
algorithms  are  complex  with  many  branches,  so  they  are  not  amenable  to  vector! zation , 
which  is  essential  for  speed  on  supercomputers.  The  pinball  algorithm  is  a 
simplified  slideline  algorithm  which  is  readily  vectorized.  Its  major  idea  is  to 
embed  pinballs  in  surface  elements  and  to  enforce  the  impenetrability  condition 
only  to  pinballs.  The  algorithm  has  been  implemented  in  the  WHAMS-3D  code. 

Examples  of  solutions  and  running  times  are  given. 


20.  DlSTRiaunON/AVAaABIUTY  OF  ABSTRACT 
□  UNCLASSIFIEOAJNUMITED  □  SAME  AS  RPT.  OOTIC  USERS 

21.  ABSTRACT  SECURITY  CLASSIFICATION 

Unclassified 

22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 

22b.  TEUFHONE  (Inchida  Atm  CodtJ  22c  OFFICE  SYMBOL 

OOfORM1473.a4MAR  «ARR*ditio«  may  until  •mn.ust.d.  SECURITY  ClASSlFiCAnON  OF  -H»S  SAGE 

Alloth«radition»;vaooio4at*.  ~ 


1 


PREFACE 

This  research  was  conducted  under  the  direction  of  co-principal  investigators 
Ted  Belytschko  and  Wing  Kam  Liu. 

The  following  students  were  supported  partially  or  in  full  by  this  project: 

Lee  P.  B  indeman 
Noreen  D.  Gilbertsen 
Edward  J.  Plaskacz 
H.Y.  Chiang 
Mark  O.  Neal 
J.S.  Chen 
Y.Y.  Lu 

The  following  degrees  were  awarded  to  graduate  students  supponed  by  the  project: 

J.  S.  Chen  (Ph.D.,  June  1989) 

Mark  O.  Neal  (Ph.D.,  December  1989) 

Noreen  D.  Gilbertsen  (Ph.D.,  December  1989) 

Edward  J.  Plaskacz  (Ph.D.,  June  1990) 

Noreen  D.  Gilbertsen  (Ph.D.,  June  1990) 

The  following  papers  were  published  as  an  outcome  of  research  supponed  by  this 

grant: 

T.  Belytschko,  X.-J.  Wang,  Z.P.  Bazant,  and  Y.  Hyun,  ‘Transient  Solutions  for  one- 
Dimensional  Problems  with  Strain  Softening,”  L  Appl.  Mech..  51(3),  513-518,  1987. 


T.J.R.  Hughes,  T.  Belytschko,  and  W.K.  Liu,  “Convergence  of  an  Element-Partitioned 
Subcycling  Algorithm  for  the  Semi-Discrete  Heat  Equation,”  Num.  Meths.  for  Partial 
Differential  Equations.  3.  131-137,  1987. 


T.  Belytschko  and  J.I.  Lin,  “A  Three-Dimensional  Impact-Penetration  Algorithm  with 


Erosion;*  Computers  &  Structures.  25(11.  95-104.  1987. 

X.-J.  Wang  and  T.  Belytschko,  “An  Efficient  Ficxurally  Superconvergent  Hexahedral 
Element,”  Engineering  Computations.  4(4),  281-288.1987. 

W.K.  Liu,  H.  Chang,  J.-S.  Chen,  and  T.  Belytschko,  “Arbitrary  Lagrangian-Eulerian 
Petrov-Galerkin  Finite  Elements  for  Nonlinear  Con tinua;*  Computer  Methods  in  Applied 
Mechanics  and  Engineering.  68.  259-310, 1988. 

J,-S.  Chen,  W.K.  Liu  and  T.  Belytschko,  “Arbitrary  Lagrangian-Eulerian  Methods  for 
Materials  with  Memory  and  friction,”  in:  Recent  Developments  in  Computational  fluid 
Dynamics.  T.  Tezduyar  and  T.J.R.  Hughes,  eds..  Applied  Mechanics  Division,  ASME, 
New  York,  AMD-Vol.  95,  1988,  11-33. 

T.  Belytschko  and  D,  Lasry,  “Nonmonotonic  Stress-Strain  Laws:  Bizarre  Behavior  and  Its 
Repercussions  on  Numerical  Solutions,”  in:  Trans.  Sixth  Army  Conf.  on  Applied 
Mathematics  and  Computing.  ARO  Report  89-1, 1989, 437-457 

T.  Belytschko,  H.S.  Chang,  and  Y.Y.  Lu,  “  A  Variationally  Coupled  Finite  Element- 
Boundary  Element  Method.**  Computers  and  Structures.  33(1).  17-20,  1989. 

T,  Belytschko  and  D.  Lasry,  “A  Study  of  Localization  Limiters  for  Strain-Softening  in 
Statics  and  Dynamics.**  Computers  and  Structures.  33(3).  707-715, 1989 

T,  Belytschko,  B.L.  Wong,  and  H.Y.  Chiang,  “Improvements  in  Low-Order  Shell 
Elements  for  Explicit  Transient  Analysis,**  in  Analytical  and  Computational  Models  of 


Shells.  A.K.  Noor,  T.  Belytschko,  and  J.C.  Simo,  eds.,  CED-Vol.  3,  ASME,  New 
York,  NY,  383-398,  1989. 


T.  Belytschko  and  M.O.  Neal,  “Contact-Impact  by  the  Pinball  Algorithm  with  Penalty, 
Projection,  and  Lagrangian  Methods,”  in  Computational  Techniques  for  Contact.  Impact. 
Penetration  and  Perforation  of  Solids.  L.E.  Schwer,  NJ.  Salamon,  and  W.K.  Liu,  cds  , 
AMD-Vol.  103,  ASME,  New  York,  NY,  97-140,  1989. 


T.  Belytschko  and  M.O.  Neal,  “The  Vectorized  Pinball  Contact  Impact  Routine,”  in  Trans. 


10th  Int.  Conf.  on  Structaral  Mechanics  in  Reactor  Technology.  A.K.  Hadjian,  ed.,Div. 


B.  lASMIRT,  Los  Angeles.  CA,  1989. 


Accession  For 
NTIS 

DTIC  TAB 
Unarm  oiniccd 
Justification- 


t 

□ 


By - 

Distribr.tir-n/ 


Avail::’-'.:!  lit  7  Codes 
|Av..::i  :..i.i/cr 
Locclai 


1 


1.  Introduction 

In  this  paper,  a  new  contact-impact  procedure  called  the  pinball  algorithm  is  described,  for 
previous  studies  of  contact-impact  see  [1-4],  A  short  description  of  pinball  was  previously  given 
by  Belytschko  and  Neal[5].  The  thrust  of  the  pinball  algorithm  is  to  allow  vectorization  of  as 
much  of  the  slideline  calculations  as  possible.  This  is  accomplished  by  greatly  simplifying  both 
the  search  for  the  elements  involved  in  the  impact  and  in  the  enforcement  of  impenetrability  witli 
the  use  of  spheres,  or  pinballs,  for  each  element  in  the  slideline  calculations.  In  this  way,  the 
search  requires  a  simple  check  on  the  distances  between  pinballs  to  determine  interpenetration. 
Once  the  contacting  pairs  of  pinballs  have  been  determined,  the  impenetrability  condition  is 
enforced  with  the  use  of  a  penalty  formulation  which  can  be  completely  vectorized.  A  similar  idea 
has  also  been  used  in  the  two-dimensional  NABOR  algorithm[6],  but  the  NABOR  method  used  an 
ad  hoc  method  based  on  spheres  for  the  determination  of  stresses  in  the  continua  and  did  not  use  a 
representation  of  the  surface  normal.  In  the  pinball  algorithm  the  element  spheres  are  used  only  in 
the  contact  algorithm,  while  standard  continuum  mechanics  is  used  for  the  continuum  elements. 

2.  Variational  Inequality  and  Discrete  Interpolants 

The  weak  form  of  the  contact  problem  is  obtained  from  the  principle  of  vinual  work  by 
appending  the  Lagrange  multiplier  type  term  5  (k  g).  We  consider  the  trial  functions  to  be 


kinematically  admissible  functions  v,  so  vj  G  V  and  A  where 

V  =  {vi:Vi6  C'’(nAu£i®),Vi  =  vf  on  i;.}  (2.1a) 

A  =  {X,:X6  C'(rc),A,<0}  (2.1b) 

As  indicated  above,  these  functions  need  only  be  piecewise  continuous  and  satisfy  essential 
boundary  conditions.  The  variations  (or  test  functions)  5vj  G  \(),  5^  G  where 

Vo  =  {  5vi :  Sv;  €  C°,  Svj  =  0  on  }  (2.2a) 

Ao  =  {5A.;5>.€  C'‘,6A.<0on  Pc}  (2.2b) 

5W  =  5W“‘  +  6M  -  5W««  (2,3) 

(see  Belytschko  [7])  where 

=  j  5v(i  dQ  (2.4) 

SI 

5we«  =  J  Svjbi  dn  +  J  SvjXj'  dr  (2.5) 
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6M  =  J  5v-pvj  dQ  (2.6) 

Q 

where  the  density  is  denoted  by  p,  the  body  force  by  bj,  the  surface  tractions  by  and  the 

Lagrange  multipliers  by  X.  The  stress  state  is  described  by  the  Cauchy  (physical)  stress  a-j. 

The  weak  form  for  the  contact  problem  is  then  given  by  : 

If  ve  V,  X  €  A  and 

5W  +  j'5{X^dr  >  0  (2.7) 

Tc 

for  all  6v  e  Vq,  6A.  €  Aq,  then  the  momentum  equation,  traction  boundary  conditions  and 
contact  surface  inequalities  are  satisfied.  The  equivalence  of  this  weak  form  to  the  governing 
partial  differential  equations  is  demonstrated  in  [5]. 

3.  Pinball  Algorithm 

The  main  idea  of  the  pinball  algorithm  is  to  enforce  the  impenetrability  condition  and 
contact  conditions  on  a  set  of  spheres,  or  pinballs  which  are  embedded  in  the  finite  elements.  By 
enforcing  the  contact  constraint  on  the  spheres  rather  than  the  elements  themselves,  the  time 
required  by  the  contact  algorithm  can  be  greatly  reduced  because:  (1)  the  determination  of  whether 
interpenetration  has  occurred  becomes  a  simple  check  of  the  distance  between  two  pinballs,  (2) 
when  combined  with  a  penalty  method,  it  involves  almost  no  recursive  calculations  or  conditional 
statements,  so  it  is  much  more  amenable  to  vectorization. 

The  hexahedral  elements  used  in  this  formulation  are  described  in  detail  by  Hanagan  and 
Belytschko[13].  A  sphere,  or  pinball,  is  embedded  in  each  of  these  hexahedral  elements  of  the 
mesh.  These  pinballs  will  then  be  used  to  determine  which  elements  are  involved  in  the  contact. 
The  center  and  radii  of  the  sphere  are  given  in  element  e  by 

Ci  =  i  i  xfj 
®  1=1 

\  47C 

respectively,  where  Cj  are  the  coordinates  of  the  center  of  the  sphere,  xj-  are  the  coordinates  of 

node  I  of  element  e,  R  is  radius  of  the  pinball,  and  V®  is  the  volume  of  element  e. 

The  center  of  each  sphere  is  simply  the  average  of  its  nodal  coordinates  while  the  radius  is 
determined  by  setting  the  volume  of  the  resulting  sphere  equal  to  the  volume  of  the  element  itself. 
For  materials  with  substantial  compressibility  the  radius  for  each  element  would  have  to  be 
recalculated  every  few  steps. 

The  detection  of  the  impacting  pairs  is,  computationally,  a  very  simple  procedure.  The 
distance  between  the  centers  of  each  slave  pinball  and  each  master  pinball  is  calculated  and  then 
compared  with  the  sum  of  the  radii  of  the  two  elements.  Interpenetration  has  occurred  when 

dij  <  Rj  +  Rj 


(3.1) 

(3.2) 


(3.3) 
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where  djj  is  the  distance  between  the  centers  of  elements  1  and  2  and  R:.  Rj  are  the  radii  of 
elements  i  and  j.  Note  that  in  this  procedure  the  masters  and  the  slaves  may  be  penetrated  by  more 
than  one  element  during  a  time  step. 

In  a  penalty  of  velocity  projection  formulation,  where  forces  are  added  to  each  node 
involved  in  contact,  contact  of  one  element  with  several  other  elements  does  not  impair  the 
algorithm.  On  the  other  hand,  in  the  displacement  projection  form  of  the  Lagrange  multiplier 
method[4],  since  the  nodes  are  moved,  difficulties  arise  when  a  single  element  comes  into  contact 
with  several  other  elements  because  nodes  are  displaced  in  the  contact  algorithm.  This  difficulty  is 
illustrated  in  figures  1.  In  figure  la  element  1  contacts  two  other  elements  2  and  3.  If  a 
displacement-based  slideline  algorithm  is  used,  element  1  would  first  be  moved  so  that  it  does  noi 
penetrate  element  2  as  in  figure  lb,  and  then  it  would  be  moved  again  as  shown  in  figure  Ic  so  as 
not  to  penetrate  element  3.  During  this  second  adjustment  element  1  may  again  penetrate  element  2 
and  this  penetration  would  go  undetected.  This  problem  is  especially  troublesome  in  vectorized 
algorithms.  Since  the  adjustment  of  position  is  a  recursive  operation,  the  first  adjustment  would 
not  take  place  at  all. 

The  penalty  force  is  proportional  to  the  depth  of  penetration  between  the  two  elements; 
therefore  the  next  step  of  the  procedure  is  to  determine  the  penetration  depth  of  the  two  elements. 
In  the  pinball  algorithm,  the  penetration  is  easily  calculated.  Consider  two  interpenetrating 
pinballs,  1  and  2,  in  figure  2,  with  the  velocities  Vj  and  V2;  the  normal  of  the  associated  surfaces 
are  and  02-  The  position  vectors  of  the  centers  of  the  two  pinballs  are  given  by  Cj  and  C2. 
The  penetration  is  given  by  g  and  is  defined  as  the  relative  displacement  of  the  centers  of  the 
pinbdls  in  the  average  normal  direction  needed  to  eliminate  interpenetration,  so  that  the  following 
equation  determines  g 

lie, -C2  +  gnl|2=  (R, +R2)2  (3.4) 

and  the  normals  of  the  two  elements  are  given  by 


n  =  j(n2-ni) 

This  gives 

g  =  -  b  +  '  c 

where 

_  ni(Xi  -  C2i) 

"j"j 

XjXj  +  C2iC2i  -  2xiC2i  -  (R,  +  Ra)^ 


(3.5) 


(5.6) 


(3.7) 


Note  that  only  the  positive  sign  on  the  radical  in  (6)  need  be  considered;  the  negative  root 
corresponds  to  a  negative  value  of  g  which  is  irrelevant. 

The  penalty  force  which  will  be  applied  to  the  nodes  of  each  element  is  proportional  to  the 
penetration  depth  and  is  given  by 

F’’  =  (Pi  g  +  P2  g)  " 


(3.8) 
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An  expression  for  the  parameter,  p2,  has  been  proposed  for  the  case  of  a  node  impacting  the 
surface  of  an  element  and  the  volume,  area  and  bulk  modulus  referred  to  the  impacted  elcment[16] 


pB 

P2  “•  V 


(3.9) 


where  B,  A,  and  V  are  the  bulk  modulus,  area  of  the  impacted  surface,  and  volume  of  the  element, 

is  the  penalty  tece  on  the  pinball. 

The  rate  of  penetration  is  computed  by 


(3.10) 


In  the  present  context,  the  properties  of  two  pinbaUs  must  be  considered,  so  the  penalty  parameter 
will  be  given  by 

P2  =  ip(BiRi+ B2R2)  (3.11) 

where  Bj,  B2  are  the  bulk  moduli  of  the  impacting  pinballs,  and  R^,  R2are  the  radii  of  the 
impacting  pinballs.  Equation  (9)  gives  the  contact  force  that  will  be  applied  in  opposite  directions 
to  each  of  the  two  impacting  pinballs.  This  force  is  then  divided  among  the  eight  nodes  of  each 
element 

=  ipP  n  =  l,8  (3.12) 

where  are  the  element  level  penalty  force  on  local  node  n  of  the  element.  These  forces  are  then 
assembled  to  the  global  force  vector  as  usual.  A  flowchart  of  the  impact  algorithm  is  given  in 
figure  3. 

The  penalty  force  is  divided  among  the  eight  nodes  of  the  hexahedron  to  preserve  the 
symmetry  of  the  underlying  linearized  system.  Since  the  position  of  the  pinball  depends  on  the 
eight  nodes  of  the  hexahedron,  the  linearized  equations  would  not  be  symmetric  if  the  force  were 
subdivided  only  among  the  surface  nodes;  an  alternative  algorithm  where  C  depends  only  on  the 
surface  node  velocities  and  hence  the  penalty  forces  are  distributed  only  to  the  surface  nodes  is 
now  under  investigation,  Belytschko  and  Bindeman[17]. 

The  penalty  forces,  along  with  the  forces  arising  from  element  stresses  and  externally 
applied  loads,  are  used  in  the  calculation  of  the  nodal  accelerations.  Therefore  the  contact  routine 
appears  in  the  algorithm  immediately  before  the  nodal  accelerations  are  calculated.  The  flowchart 
of  the  complete  explicit  algorithm  with  the  contact  algorithm  are  given  in  figure  4. 

4.  Numerical  Examples 

In  order  to  test  the  accuracy  and  efficiency  of  the  proposed  contact  procedure  several 
example  problems  were  examined.  The  first  problem  considered  was  the  impact  of  two  one¬ 
dimensional  bars  consisting  of  ten  elements  each.  This  problem  was  performed  in  order  to 
compare  two  different  methods  of  enforcing  the  impenetrability  constraint:  the  penalty  method  and 
the  projection  method.  This  contact  constraint  is  the  only  nonlinearity  that  appears  in  the  problem. 
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As  can  be  seen  in  figure  5  one  of  the  bars  is  given  an  initial  velocity  so  that  it  impacts  with  the 
second  bar.  The  material  properties  are  such  that  the  wave  speed  in  the  two  bars  is  10.0  m/s. 

Figure  6  through  8  gi’^e  the  velocity  time  histories  for  the  nodes  at  the  midpoint  of  the  first 
rod  (x=5.0),  at  the  interface  on  the  first  rod  (x=10.0),  and  at  the  midpoint  of  the  second  rod 
(x=15.5).  As  can  be  seen  from  figure  7,  the  penalty  method  gives  a  rather  noisy  solution  at  the 
contact  interface  yet  this  does  not  appear  to  have  much  affect  away  from  the  contact  zone  (see 
figures  6  and  8).  The  results  for  the  Lagrange  multiplier  method  were  nearly  identical  to  the 
projection  metiiod  and  therefore  were  not  included  in  the  results. 

The  second  problem  considered  was  of  a  copper  rod  impacting  a  steel  plate  at  high 
velocity.  The  rod,  acting  as  the  projectile,  consisted  of  414  elements  while  the  plate  or  target  was 
made  up  of  1,428  elements.  The  geometry  and  the  material  properties  for  each  of  the  objects  is 
given  in  table  1.  The  evolution  of  the  problem  is  shown  in  figure  9.  In  this  example  problem  and 
the  one  that  follows  a  Von  Mises  yield  criteria  is  used  with  a  piece- wise  linear  stress-strain  model. 
Once  the  effective  stress  reaches  yield  the  material  undergoes  isotropic  hardening  until  the  effective 
stress  is  equal  to  the  ultimate  stress,  at  which  point  the  material  is  considered  perfectly  plastic. 
When  the  effective  strain  of  an  element  reaches  the  maximum  allowable  effective  strain,  the 
element  is  eroded,  that  is  the  stress  in  the  element  is  considered  zero  fi'om  that  time  on.  The 
maximum  allowable  effective  strain  used  for  steel  is  1.0  while  that  used  for  copper  is  2.0. 

This  problem  was  also  examined  by  Belytschko  and  Lin[4]  with  their  projection  method 
and  a  comparison  of  running  times  for  both  the  methods  is  given  in  the  first  column  of  table  2.  As 
can  be  seen  firom  this  table,  the  efficiency  of  the  new  algorithm  is  substantially  better  than  the 
previous  one  on  a  vectorized  machine.  For  this  comparison,  both  algorithms  were  implemented 
into  the  three  dimensional  finite  element  code  WHAM3D  and  run  on  a  Cray  X-MP/14  with  the 
CFT77  compiler.  The  differences  shown  in  running  times  is  due  only  to  the  different  slideline 
algorithms  used. 

The  importance  of  vectoiization  for  both  of  these  slideline  algorithms  is  demonstrated 
graphically  in  figures  10  and  1 L  In  these  figures  a  breakdown  of  the  total  CPU  is  given  for  each 
part  of  the  program  for  both  vectorized  and  unvectorized  compilers.  (The  unvectorized  runs  were 
performed  on  tiie  same  machine  using  the  CFT77  compiler  with  the  vectoiization  turned  off.)  For 
unvectorized  mns,  the  new  algorithm  is  only  marginally  more  efficient  than  the  previous  method. 
When  the  vectorized  compiler  is  used,  however,  the  old  version  of  the  slideline  algorithm 
consumes  nearly  fifty  percent  of  the  total  CPU,  the  new  procedure  this  value  has  been  reduced  to 
only  fifteen  percent  for  the  vectorized  run. 

The  impact  of  an  elastic  sphere  with  a  rigid  wall  was  examined  to  test  the  accuracy  of  the 
new  algorithm.  Ten  degrees  of  the  sphere  are  modelled  with  499  elements  as  shown  in  figure  12. 
All  of  the  nodes  in  this  model  are  constrained  in  the  circumferential  direction;  in  addition,  the 
nodes  along  the  diameter  are  also  constrained  in  the  radial  direction.  Figure  12  also  gives  the 
material  properties  and  dimensions  of  the  sphere.  The  contact  radius  as  a  function  of  time  is 
compared  for  the  numerical  simulation  and  the  analytical  result[21]  in  figure  13. 

The  final  example  problem  is  that  a  rod  striking  a  plate  at  136,000  cm/sec  and  a  60^ 
obliquity.  The  projectile  is  modeled  with  414  elements  and  the  target  with  5,300  elements.  The 

time  duration  of  the  simulation  is  1,0  x  lO’*^  sec.  The  geometry  and  material  properties  are 
described  in  table  3.  The  initial  mesh  is  shown  in  figure  14.  The  final  mass  of  the  projectile  is 
42%  of  the  initial  mass  and  the  exit  velocity  is  1 17,756.34  cm/sec  in  the  x  direction. 

5.  Conclusions 

The  major  breakthrouth  of  this  paper  is  the  demonstration  that  a  contact-impact  algorithm 
can  be  simplified  dramatically  by  interpreting  the  gap  g  between  the  bodies  as  the  gap  between 
spheres  embedded  in  the  elements.  This  simplifies  the  contact-impact  algorithm  and  facilitates 
vectoiization.  Computer  times  for  large  three-dimensional  problems  show  a  fivefold  speedup  in 
the  slideline  algorithm  and  as  much  as  a  factor  of  two  in  total  running  time 
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Table  1.  Geometry  and  material  properties  of  penetration  problem  1. 


ProjcgtilcCEod  with  a  round  iipsg) 

Tareeti'plate') 

Length 

= 

4.900  in. 

3.950  in. 

Width 

= 

- 

7.900  in.fhalf  plate  is  modeled) 

Thickness 

= 

- 

0.375  in. 

Radius 

= 

0.500  in. 

- 

Density 

= 

8.31e-3  lb-sec^  /  in.^ 

7.34e-3  lb-sec2  /  in .4 

Bulk  Modulus 

= 

2.0739e+7  psi 

2.4200e+7  psi 

Shear  Modulus 

= 

6.3800e+6  psi 

9.3000e+6  psi 

Plastic  Modulus 

= 

1.5000e+5  psi 

1.4300e+5  psi 

Yield  Stress 

= 

2.0300e+4  psi 

1.6000e+5  psi 

Ultimate  Stress 

= 

6.5300e*f4  psi 

2.0300e+5  psi 

Initial  Velocity 

= 

5.5566e+4  in./  sec.(x-component) 

0.0 

-5.5566e+4  in./  sec  (z-component) 


Table  2.  Timing  studies  for  penetration  problems. 

Algorithm  Example  1  Example  2.  mesh  1  Example  2.  mesh  2 

Previous  method  34,7  sec,  94.4  sec.  302.2  sec. 

Pinball  algorithm  22.0  sec.  '0.4  sec.  143.0  sec. 


Table  3.  Geometry  and  material  properties  of  penetration  problem  4 


Length 

— 

Proiectilefrod  with  a  round  nose) 

10.25  cm 

Tareetfulate) 

20.00cm 

With 

=. 

- 

10.00cm 

Thickness 

=. 

- 

2.53cm 

Radius 

= 

0.51cm 

- 

Density 

7.77e-3  kgm/cm^ 

7.77e-3  kgm/cm^ 

Young’s  Modulus 

2.07e+7  N/cm2 

2.07e+7  N/cm2 

Yield  Stress 

= 

1.38e+5  N/cm2 

9.13e+4  N/cm2 

Ultimate  stress 

- 

1.12e+5  N/cm2 

Failure  stress 

= 

- 

3.65e+4  N/cm2 

Failure  Strain 

= 

2.0 

2.5 

Initial  Veloaty 

1 17756.34  cm/sec.(  x-component) 
-68015.79cm/sec  ( y-component) 

Fig,  2.  Interpenetration  of 
two  pinballs. 


Fig.  1.  a)  Impact  of  one  element  with  two  other  elements; 

b)  movement  of  element  1  due  to  impact  with  element  2; 

c)  movement  of  element  1  due  to  impact  with  clement  3. 
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1 .  If  this  is  the  first  step,  use  the  element  volume  to  calculate  a  radius  for  all 
elements  on  the  slideline. 

2.  Calculate  element  normals  for  all  elements.  Elements  with  zero  normals  are 
eliminated  from  consideration  in  the  contaa  search. 

3 .  Calculate  the  center  of  all  elements  with  non-zero  normals. 

4 .  Put  elements  into  appropriate  cells. 

5 .  Loop  through  elements  of  each  cell  to  determine  the  penetrating  pairs  of 
elements. 

6 .  Calculate  the  contact  forces  to  be  applied  to  the  nodes  of  impacting  element 
pairs. 

7 .  Return  to  main  driving  routine. _ 


Fig.  3.  Pinball  algorithm. 


1 .  Initialization. 

2 .  Calculate  the  external  nodal  forces. 

3.  Compute  the  internal  nodal  force  array. 

a.  Calculate  the  element  stresses. 

b.  Compute  the  element  nodal  forces  arising  from  the  element  stresses. 

c.  Assemble  the  element  nodal  forces  to  the  internal  nodal  force  array. 

4.  Call  the  slidehne  algorithm  to  calculate  the  contact  forces  and  add  them  to  the 
external  force  array. 

5.  Compute  the  nodal  accelerations. 

6.  Integrate  the  accelerations  to  obtain  the  nodal  velocities  and  displacements. 

7.  Go  to  2. 


Fig.  4.  Explicit  time  stepping  procedure  including  slideline  procedure. 


velocity 


Fig.  5.  One-dimensional  impact  problem. 


0.0  3.0  6.0  9.0  12.0 

Time(sec.) 


Fig.  6.  Projection  and  penalty  methods  at  midpoint  of  first  rod. 


Velocily 


0.0 


3.0 


9.0 


12.0 


6.0 

Time(sec.) 

Fig.  7.  Projection  and  penalty  methods  at  interface  on  first  rod. 
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Fig.  8.  Projection  and  penalty  methods  at  midpoint  of  second  rod. 


1  2 


Fig.  9.  Example  problem  1  at  times  0,  25,  and  54  pseconds. 


Fig.  10.  CPU  rcquiremenis  of  the  Belytschko-Lin  algorithm 


CPU  CPU 


1  4 


Fig.  11.  CPU  requirements  of  the  pinball  algorithm. 


Contact  Radius  (inches) 


E  =  ICXX)  psi  p  =  0.0 1  ibf-secVin'* 

\)  =  0.3  Vq  =  2.0  in/scc. 

Radius  =  5.0  in 

Fig.  12.  Elastic  impact  of  sphere  with  rigid  wall. 


Time  (seconds) 

Fig.  13.  Radius  of  contact. 


Fig.  14.  The  initial  mesh  of  a  rod  striking  a  plate. 
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